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ABSTRACT 

We present a new method for calculating arrival distribution of Ultra-High Energy Cosmic Rays (UHECRs) 
including modifications by the galactic magnetic field. We perform numerical simulations of UHE anti-protons, 
which are injected isotropically at the earth, in the Galaxy and record the directions of velocities at the earth and 
outside the Galaxy for all of the trajectories. We then select some of them so that the resultant mapping of the 
velocity directions outside the Galaxy of the selected trajectories corresponds to a given source location scenario, 
applying Liouville's theorem. We also consider energy loss processes of UHE protons in the intergalactic space. 
Applying this method to our source location scenario which is adopted in our recent study and can explain the 
AGASA observation above 4 x 10 19 eV, we calculate the arrival distribution of UHECRs including lower energy 
(E > 10 19 eV) ones. We find that our source model can reproduce the large-scale isotropy and the small-scale 
anisotropy on UHECR arrival distribution above 10 19 eV observed by the AGASA. We also demonstrate the 
UHECR arrival distribution above 10 19 eV with the event number expected by future experiments in the next 
few years. The interesting feature of the resultant arrival distribution is the arrangement of the clustered events 
in the order of their energies, reflecting the directions of the galactic magnetic field. This is also pointed out by 
Alvarez-Muniz, Engel & Stanev (2002). This feature will allow us to obtain some kind of information about the 
composition of UHECRs and the magnetic field with increasing amount of data. 

Subject headings: cosmic rays — methods: numerical — ISM: magnetic fields — galaxies: general — 
large-scale structure of universe 



1. INTRODUCTION 

There is no statistically significant large scale anisotropy in 
the observed arrival distribution of ultra-high energy cosmic 
rays (UHECRs) above 10 19 eV (Takeda et al. 1999). This may 
imply an extragalactic origin of cosmic rays above 10 19 eV, 
combined with the change of spectral slope of the observed 
energy spectrum at ~ 10 19 eV (Bird et al. 1994; Yoshida et 
al. 1995; Takeda et al. 1998). Another important feature of 
the UHECR arrival distribution is the small scale clusterings 
of the arrival directions (Takeda et al. 1999, 2001). The cur- 
rent AGASA data set of 57 events above 4 x 10 19 eV contains 
four doublets and one triplet within a separation angle of 2.5°. 
Chance probability to observe such clusters under an isotropic 
distribution is only about 1 % (Hayashida et al. 2000; Takeda et 
al. 2001). 

On the other hand, the cosmic -ray energy spectrum does 
not show the GZK cutoff (Greisen 1966; Zatsepin & Kuz'min 
1966) because of photopion production with the photons of the 
cosmic microwave background (CMB) and extends above 10 20 
eV (Takeda et al. 1998). The discrepancy between the AGASA 
and the High Resolution Fly's Eye (HiRes; Wilkinson et al. 
1999), which reports the cosmic ray flux with the GZK cut-off 
around 10 20 eV (Abu-Zayyad et al. 2002), remains to be one 
of the major open question in astroparticle physics. This is- 
sue is left for future investigation by new large-aperture detec- 
tors under development, such as South and North Auger project 
(Capelle et al. 1998), the EUSO (Benson & Linsley 1982), and 
the OWL (Cline & Stecker 2000) experiments. 

In our recent work (Yoshiguchi et al. 2003a, hereafter Paper 
I), we perform numerical simulations for propagation of UHE 
protons in intergalactic space, and examine whether the present 



AGASA observation above 4 x 10 19 eV can be explained by 
a bottom-up scenario in which the source distribution of UHE- 
CRs is proportional to that of galaxies. We use the Optical Red- 
shift Survey (ORS; Santiago et al. 1995) to construct realistic 
source models of UHECRs. 

In Paper I, we calculate both the energy spectrum and arrival 
directions of UHE protons, and compare the results with the 
AGASA observation above 4 x 10 19 eV. We find that the large- 
scale isotropy and the small-scale anisotropy of the UHECR 
arrival distribution observed by the AGASA can be reproduced 
when ~ 1/50 of the ORS sample more luminous than -20.5 
mag are selected as UHECR sources, in the case of weak ex- 
tragalactic magnetic field (EGMF B < 1 nG). In terms of the 
source number density, this constraint corresponds to ~ 10~ 6 
Mpc" 3 . 

The small scale anisotropy can not be well reproduced in 
the case of strong EGMF (B > 10 nG), because the correla- 
tion at small scale between events which originate from a sin- 
gle source is eliminated, or the correlation continues to larger 
angle scale, due to large deflection when UHECRs propagate 
in the EGMF from sources to the earth. Although Isola & Sigl 
(2002) and Sigl, Miniati, & Ensslin (2003) conclude that the 
expected small-scale anisotropy and large-scale isotropy for lo- 
cal enhancement of UHECR sources in the LSC in the presence 
of the strong EGMF (~ 1/i G) are in marginal agreement with 
the AGASA, the consistency is somewhat worse than that pre- 
dicted by our scenario for B = 1 nG. Of course, we can not draw 
any firm conclusion about the strength of the EGMF, consider- 
ing the current limited amount of data. However, we assume 
extremely weak EGMF throughout the paper. 

If local enhancement of UHECR sources in the LSC (Isola 
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& Sigl 2002; Sigl, Miniati, & Ensslin 2003) is disfavored from 
the observations, there is no way that explains the observed 
extension of the cosmic-ray spectrum beyond the GZK cutoff. 
Our conclusion in Paper I is that a large fraction of cosmic rays 
above 10 20 eV observed by the AGASA experiment might orig- 
inate in the top-down scenarios, or that the energy spectrum 
measured by the Hires experiment might be better. 

As mentioned above, we obtain the constraint on the source 
number density as ~ 10~ 6 Mpc~ 3 by comparing our model pre- 
diction with the AGASA data only above 4 x 10 19 eV. It is very 
important to examine whether our source model can explain the 
AGASA data including lower energy (~ 10 19 eV) one. On the 
other hand, the arrival directions of UHECRs above 10 19 eV are 
modified by the galactic magnetic field (GMF) by a few - ~ 10 
degrees. In order to accurately calculate the expected UHECR 
arrival distribution and compare with the observations, the ef- 
fect of the GMF should be taken into account. 

The first step of the studies of UHECR propagation in the 
GMF is found in Alvarez-Muniz, Engel & Stanev (2002). 
Alvarez-Muniz, Engel & Stanev (2002) calculate the expected 
arrival distribution of UHECRs above 10 19 4 eV for several 
source location scenarios. They perform numerical simulations 
of UHECR propagation in the Galaxy injected from sources 
toward the earth. The radius of the earth (detector) must be 
so small that the unavoidable smearing in arrival angle is kept 
smaller than the accuracy of arrival direction determination ~ a 
few degree (Takeda et al. 2001). In this case, the number frac- 
tion of injected UHECRs arriving at the earth is very small. 
This requires a large number of particles to be propagated, 
which takes enormous CPU time. 

In this paper, we present a new method for calculating 
UHECR arrival distribution which can be applied to several 
source location scenarios including modifications by the GMF. 
We numerically calculate the propagation of anti-protons from 
the earth toward the outside of the Galaxy (in this study, we 
set a sphere centered around the Galactic center with radius 
?"src =40 kpc as the boundary condition), including the effects 
of Lorentz force due to the GMF. The anti-protons are ejected 
isotropically from the earth. By this calculation, we can obtain 
the trajectories and the sky map of position of anti-protons that 
have reached to the boundary at radius r src = 40 kpc. 

Next, we regard the obtained trajectories as the ones of PRO- 
TONS from the outside of the galaxy toward the earth. Also, 
we regard the obtained sky map of position of anti-proton at 
the boundary as relative probability distribution (per steradian) 
for PROTONs to be able to reach to the earth for the case in 
which the flux of the UHE protons from the extra-galactic re- 
gion is isotropic (in this study, this flux corresponds to the one 
at the boundary r src = 40 kpc). This treatment is supported by 
the Liouville's theorem. When the flux of the UHE protons at 
the boundary is anisotropic (e.g., the source distribution is not 
isotropic), this effect can be included by multiplying this ef- 
fect (that is, by multiplying the probability density of arrival 
direction of UHE protons from the extra-galactic region at the 
boundary) to the obtained relative probability density distribu- 
tion mentioned above. 

By adopting this new method, we can consider only the tra- 
jectories of protons which arrive to the earth, which, of course, 
helps us to save the CPU time efficiently and makes calcula- 
tion of propagation of CRs even with low energies (~ 10 19 eV) 
possible within a reasonable time. We also consider the energy 
loss processes of UHE protons in the intergalactic space, which 



is not taken into account by Alvarez-Muniz, Engel & Stanev 
(2002). 

With this method, we calculate the UHECR arrival distribu- 
tion above 10 19 eV for our source scenario which can explain 
the current AGASA observation above 10 19 6 eV. Using the har- 
monic amplitude and the two point correlation function as sta- 
tistical quantities, we compare our model prediction with the 
AGASA observation. We also demonstrate the arrival distri- 
bution of UHECRs with the event number expected by future 
experiments such as South and North Auger project (Capelle et 
al. 1998), the EUSO (Benson & Linsley 1982), and the OWL 
(Cline & Stecker 2000) experiments. 

In section 2, we introduce the GMF model. We explain the 
method for calculating UHECR arrival distribution in section 3. 
Results are shown in section 4. In section 5, we summarize the 
main results. 

2. GALACTIC MAGNETIC FIELD 

In this study, we adopt the GMF model used by Alvarez- 
Muniz, Engel & Stanev (2002), which is composed of the spiral 
and the dipole field. In the following, we briefly introduce this 
GMF model. 

Faraday rotation measurements indicate that the GMF in the 
disk of the Galaxy has a spiral structure with field reversals at 
the optical Galactic arms (Beck 2001). We use a bisymmetric 
spiral field (BSS) model, which is favored from recent work 
(Han, Manchester, & Qiao 1999; Han 2001). The Solar Sys- 
tem is located at a distance n = /?® =8.5 kpc from the center 
of the Galaxy in the Galactic plane. The local regular mag- 
netic field in the vicinity of the Solar System is assumed to be 
fisoiar ~ 1.5 jiG in the direction I = 90°+ p where the pitch angle 
is p = -10° (Han & Qiao 1994). 

In the polar coordinates (r\ \ , <f)), the strength of the spiral field 
in the Galactic plane is given by 



B(r||»= J B ^-f-J cos^-/31n-^j (1) 

where B Q = 4.4 fj,G, r Q = 10.55 kpc and (3 = 1 /tmp = -5.67. 
The field decreases with Galactocentric distance as \/r\\ and it 
is zero for r\\ > 20 kpc. In the region around the Galactic center 
(ni < 4 kpc) the field is highly uncertain, and thus assumed to 
be constant and equal to its value at r| | =4 kpc. 

The spiral field strengths above and below the Galactic plane 
are taken to decrease exponentially with two scale heights 
(Stanev 1996) 

ibc a m \j>t ,am / ex P(" z ) : |z|<0.5kpc 
\B{r\\A,z)\ = \B{r^)\[ exp( z3) exp(?) : | z | > 0.5 kpc 

(2) 

where the factor exp(-3/8) makes the field continuous in z. The 
BSS spiral field we use is of even parity, that is, the field direc- 
tion is preserved at disk crossing. 

Observations show that the field in the Galactic halo is much 
weaker than that in the disk. In this work we assume that the 
regular field corresponds to a AO dipole field as suggested in 
(Han 2002). In spherical coordinates (r,6,(p) the (x,y,z) com- 
ponents of the halo field are given by: 

B x = -3 /ig sin#cos#cos(/3/r 3 

By = —3 /iG sin cos 9 sin ip/r 3 (3) 

B z = fiQ (1- 3 cos 2 9)/r 3 

where ^g ~ 184.2 ^G kpc 3 is the magnetic moment of the 
Galactic dipole. The dipole field is very strong in the central 
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region of the Galaxy, but is only 0.3 pG in the vicinity of the 
Solar system, directed toward the North Galactic Pole. 

There is a significant turbulent component, B ran , of the Galac- 
tic magnetic field. Its field strength is difficult to measure and 
results found in literature are in the range of B mn = 0.5 . . . 2B Kg 
(Beck 2001). However, we neglect the random field throughout 
the paper, in order to make easy to see the effects of the regular 
field, such as the arrangement of the clustered event in the or- 
der of their energies (section 4.3). Possible dependence of the 
results on the random field is discussed in the section 4.2. 



3. NUMERICAL METHOD 
3.1. Propagation of UHECRs in the Intergalactic Space 

The energy spectrum of UHECRs injected at extragalactic 
sources is modified by the energy loss processes when they 
propagate in the intergalactic space. This subsection provides 
the method of Monte Carlo simulations for propagation of UHE 
protons in intergalactic space. 

We assume that UHECRs are protons injected with a power 
law spectrum within the range of (10 19 - 10 22 )eV. 10000 pro- 
tons are injected in each of 31 energy bins, that is, 10 bins per 
decade of energy. Then, UHE protons are propagated includ- 
ing the energy loss processes (explained below) over 3 Gpc for 
15 Gyr. We take a power law index as 2.6 in order to fit the 
calculated energy spectrum to the one observed by the AGASA 
(Marco, Blasi, & Olinto 2003). 

UHE protons below ~ 8 x 10 19 eV lose their energies mainly 
by pair creations and adiabatic energy losses, and above it by 
photopion production (Berezinsky & Grigorieva 1988; Yoshida 
& Teshima 1993) in collisions with photons of the CMB. We 
treat the adiabatic loss as a continuous loss process. We cal- 
culate the redshift z of source at a given distance using the 
cosmological parameters Ho = 71 km s" 1 Mpc -1 , £l m = 0.27, 
and £1a = 0.73. Similarly, the pair production can be treated 
as a continuous loss process considering its small inelasticity 
(~ 10~ 3 ). We adopt the analytical fit functions given by Chodor- 
owski, Zdziarske, & Sikora (1992) to calculate the energy loss 
rate for the pair production on isotropic photons. The same ap- 
proach has been adopted in our previous studies (Paper I, Ide 
et al. 2001 ; Yoshiguchi, Nagataki, & Sato 2003c). 

On the other hand, protons lose a large fraction of their en- 
ergy in the photopion production. For this reason, its treatment 
is very important. We use the interaction length and the en- 
ergy distribution of final protons as a function of initial proton 
energy which is calculated by simulating the photopion produc- 
tion with the event generator SOPHIA (Mucke et al. 2000). 

In this study, we neglect the effect of the EGMF because 
of the following two reasons. First, numerical simulations 
of UHECR propagation in the EGMF including lower energy 
(~ 10 19 eV) ones take a long CPU time. Second, we show in 
our previous study that small scale clustering can be well re- 
produced in the case of weak EGMF (B < InG) (Paper I). Isola 
& Sigl (2002) and Sigl, Miniati, & Ensslin (2003) show that 
the expected small scale anisotropy for local enhancement sce- 
nario of UHECR sources in the presence of strong EGMF (~ 1 p 
G) in the Local Super Cluster is marginally consistent with the 
AGASA observation. However, the consistency of small-scale 
anisotropy is somewhat worse than that predicted by our sce- 
nario in the case of weak EGMF (Paper I). Although we can 
not draw any firm conclusion because of the limited amount of 
data, we assume extremely weak EGMF throughout the paper. 



3.2. Source Distribution 

In this study, we apply the method for calculating the 
UHECR arrival distribution with modifications by the GMF 
(section 3.3) to our source location scenario, which is con- 
structed by using the ORS (Santiago et al. 1995) galaxy catalog. 
As mentioned in section 1 , we show in Paper I that the arrival 
distribution of UHECRs observed by the AGASA can be repro- 
duced when ~ 1 /50 of the ORS galaxies more luminous than 
M\ lm = -20.5 is selected as UHECR sources. We consider only 
this source model throughout the paper. It is unknown how 
much an ultimate UHECR source contribute to the observed 
cosmic ray flux. In paper I, we thus consider the two cases in 
which (1) all galaxies inject the same amount of cosmic rays, or 
(2) they inject cosmic rays proportionally to their absolute lu- 
minosity. However, we find that the results in the two cases do 
not differ from each other, as far as we focus on the luminous 
galaxies as UHECR sources. Accordingly, we restrict ourselves 
to the case that all galaxies inject the same amount of cosmic 
rays. 

In order to calculate the energy spectrum and the distribu- 
tion of arrival directions of UHECRs realistically, there are two 
key elements of the galaxy sample to be corrected. First, galax- 
ies in a given magnitude-limited sample are biased tracers of 
matter distribution because of the flux limit (Yoshiguchi et al. 
2003b). Although the sample of galaxies more luminous than 
-20.5 mag is complete within 80 h~ l Mpc (where h is the Hub- 
ble constant divided by 100 km s" 1 and we use h = 0.71), it 
does not contain galaxies outside it for the reason of the se- 
lection effect. We distribute sources of UHECRs outside 80 h~ l 
Mpc homogeneously. Their number density is set to be equal to 
that inside 80 h~ l Mpc. We do not take into account luminosity 
evolution for simplicity. 

Second, our ORS sample does not include galaxies in the 
zone of avoidance {\b\ < 20°). In the same way, we distribute 
UHECR sources in this region homogeneously, and calculate 
its number density from the number of galaxies in the observed 
region. 

3.3. Calculation of the UHECR Arrival Distribution with 
modifications by the GMF 

In this subsection, we present the method of calculation of 
UHECR arrival distribution with modifications by the GMF. We 
start by injecting anti-protons from the earth isotropically, and 
follow each trajectory until 

1 . anti -proton reaches a sphere of radius 40 kpc centered at 
the galactic center, or 

2. the total path length traveled by anti-proton is larger than 
200 kpc, 

by integrating the equations of motion in the magnetic field. 
It is noted that we regard these anti-protons as PROTONs in- 
jected from the outside of the Galaxy toward the earth. The 
number of propagated anti-proton is 2,000,000. We have 
checked that the number of trajectories which are stopped by 
the limit (2) is smaller than 0.1% of the total number. The en- 
ergy loss of protons can be neglected for these distances. Ac- 
cordingly, we inject anti-protons with injection spectrum simi- 
lar to the observed one ~ E~ 2 1 . (Note that this is not the energy 
spectrum injected at extragalactic sources.) 

The result of the velocity directions of anti-protons at the 
sphere of radius 40 kpc is shown in the right panel of figure 1 
in the galactic coordinate. From Liouville's theorem, if the 
cosmic-ray flux outside the Galaxy is isotropic, one expects an 
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isotropic flux at the earth even in the presence of the GMF. This 
theorem is confirmed by numerical calculations shown in figure 
6 of Alvarez-Muniz, Engel & Stanev (2002), which is the same 
figure as our figure 1 except for threshold energy. Thus, the 
mapping of the velocity directions in the right panel of figure 1 
corresponds to the sources which actually give rise to the flux 
at the earth in the case that the sources (including ones which 
do not actually give rise to the flux at the earth) are distributed 
uniformly. 

We calculate the UHECR arrival distribution for our source 
scenario using the numerical data of the propagation of UHE 
anti-protons in the Galaxy. Detailed method is as follows. At 
first, we divide the sky into a number of bins with the same 
solid angle. The number of bins is taken to be 360(0 x 200(b). 
We then distribute all the trajectories into each bin according 
to their directions of velocities (source directions) at the sphere 
of radius 40 kpc. Finally, we randomly select trajectories from 
each bin with probability P se iec(j,k, E) defined as 



1 dN/dE(di,E) 
df ' 



(4) 



Here subscripts j and k distinguish each cell of the sky, di is 
distance of each galaxy within the cell of (j,k), and the summa- 
tion runs over all of the galaxies within it. E is the energy of 
proton, and dN/dE(dj,E) is the energy spectrum of protons at 
our galaxy injected at a source of distance di. 

The normalization of P se i sc (j,k,E) is determined so as to set 
the total number of events equal to a given number, for example, 
the event number of the current AGASA data. When f se i ec > 1, 
we newly generate events with number of (P se iec- 1) x N(j,k), 
where N(j,k) is the number of trajectories within the sky cell 
of (j,k). The arrival angle of newly generated proton (equiva- 
lently, injection angle of anti-proton) at the earth is calculated 
by adding a normally distributed deviate with zero mean and 
variance equal to the experimental resolution 2.8° (1.8°) for 
E > 10 19 eV (4 x 10 19 eV) to the original arrival angle. We per- 
form this event generation 20 times in order to calculate the av- 
erages and variances, due to the finite number of the simulated 
events, of the statistical quantities (section 3.4). 

3.4. Statistical Methods 

In this subsection, we explain the two statistical quantities, 
the harmonics analysis for large scale anisotropy (Hayashida 
et al. 1999), the two point correlation function for small scale 
anisotropy. 

The harmonic analysis to the right ascension distribution 
of events is the conventional method to search for global 
anisotropy of cosmic ray arrival distribution. For a ground- 
based detector like the AGASA, the almost uniform observation 
in right ascension is expected. The m-th harmonic amplitude r 
is determined by fitting the distribution to a sine wave with pe- 
riod 2ir/m. For a sample of n measurements of phase, fa, fa, 
■■-,(/)„ (0 < <j>i < 27r), it is expressed as 

r = (a 2 + b 2 ) 1/2 (5) 

where, a = -S" =1 cosmfa, b = -S" =1 sinm^. We calculate the 
harmonic amplitude for m = 1 , 2 from a set of events generated 
by the method explained in the section 3.3. 

If events with total number n are uniformly distributed in 
right ascension, the chance probability of observing the ampli- 
tude > r is given by, 

P = exp(-fe), (6) 



where 

k = nr 1 /4. (7) 

The current AGASA 775 events above 10 19 eV is consistent 
with isotropic source distribution within 90 % confidence level 
(Takeda et al. 2001). We therefore compare the harmonic am- 
plitude for P = 0.1 with the model prediction. 

The two point correlation function N(8) contains information 
on the small scale anisotropy. We start from a set of generated 
events. For each event, we divide the sphere into concentric 
bins of angular size A8, and count the number of events falling 
into each bin. We then divide it by the solid angle of the corre- 
sponding bin, that is, 



N(0) = 



E 



1 [sr- 1 ], (8) 



27r|cos6»-cos((9+A(9)| 

9<<p<9+A9 

where <f> denotes the separation angle of the two events. AO 
is taken to be 1° in this analysis. The AGASA data shows 
correlation at small angle (~ 3°) with 2.3 (4.6) a significance 
of deviation from an isotropic distribution for E > 10 19 eV 
(£ > 4 x 10 19 eV) (Takeda et al. 2001). 

4. RESULTS 

4.1. Arrival Distribution ofUHECRs above 10 19 eV 

In this subsection, we present the results of the arrival distri- 
bution of UHECRs above 10 19 eV, using the method explained 
in the section 3.3. At first, figure 2 shows the distribution of the 
sources for a specific source selection when ~ 1 /50 of the ORS 
galaxies more luminous than Mr, m = -20.5 is randomly selected 
as UHECR sources in the galactic coordinate. We show only 
the sources within 300 Mpc from us for clarity, as circles of 
radius inversely proportional to their distances. It is noted that 
the sources outside 1 13(= 80/z -1 ) Mpc are randomly distributed 
because the ORS sample does not contain any galaxy outside it. 
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FIG. 3. — Energy spectrum with injection spectrum E~ 26 , predicted by 
the source model of figure 2. The contributions from sources at different dis- 
tances are also shown. We also show the observed cosmic-ray spectrum by the 
AGASA experiments (Hayashida et al. 2000). 

We show in figure 3 the expected energy spectrum for the 
source model of figure 2. The injection spectrum is set to be 
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FIG. 1. — Arrival directions of anti-protons with E > 10 19 eV at the sphere of Galactocentric radius of 40 kpc (right panel) in the galactic coordinate, 
anti-protons are injected at the earth isotropically (as shown in the left panel) with an injection spectrum E~ 2J . 
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FIG. 2. — Distribution of the sources in our model in the galactic coordinate. We show only the sources within 300 Mpc from us as circles of radius inversely 
proportional to their distances. It is noted that the sources outside 113(= 80/T 1 ) Mpc are randomly distributed because the ORS sample does not contain any galaxy 
outside it. 
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E~ 2 6 . The contributions from sources at different distances are 
also shown. We also show the observed cosmic -ray spectrum 
by the AGASA experiments (Hayashida et al. 2000). The re- 
sultant spectrum is in good agreement with the one observed by 
the AGASA, except for E > 10 20 eV. As mentioned above, we 
conclude in Paper I that a large fraction of cosmic rays above 
10 20 eV observed by the AGASA experiment might originate in 
the top-down scenarios. Accordingly, we consider only cosmic 
rays with E < 10 20 eV throughout the paper. 

As mentioned above, Alvarez-Muniz, Engel & Stanev (2002) 
does not take the energy loss processes in the intergalactic space 
into account. Thus, they can not include the effects of dif- 
ference between resultant energy spectra injected at different 
distances into numerical calculations. In our calculations, how- 
ever, sources at larger distance mainly contribute to the cosmic 
ray flux at lower energies as is evident from figure 3. This en- 
ables us to calculate the arrival distribution of UHECRs under 
more realistic situations. 

Given the source distribution and the resultant energy spec- 
trum as a function of the source distance, we can calculate the 
right hand side of Eq. 4. Then we perform the selection of tra- 
jectories according to the probability fseiec, as explained in the 
section 3.3. 

One realization of the event generations is shown in figure 4. 
The events are shown by color according to their energies. This 
figure corresponds to figure 1. That is, the injection directions 
of anti -protons at the earth (figure 1, left) corresponds to the ar- 
rival directions of protons (figure 4, left). Similarly, the arrival 
directions of anti-protons at the sphere of Galactocentric radius 
of 40 kpc (figure 1, right) does to the directions of the sources 
which actually give rise to the cosmic -ray flux (figure 4, right). 

For the source model of figure 2, the nearest source is lo- 
cated at (b,l) = (31°, 284°) and 64 Mpc from us. A number of 
the simulated events are clustered at this direction as seen in 
figure 4. Furthermore, these events are aligned in the sky ac- 
cording to the order of their energies, reflecting the direction of 
the GMF at this direction. As we will show in the section 4.3, 
this interesting feature of the UHECR arrival distribution be- 
comes evident with increasing amount of the event number. 

4.2. Statistics on the UHECR Arrival Distribution 

In this subsection, we show the results of the statistical quan- 
tities on the UHECR arrival distribution above 10 19 eV. In the 
last section, we showed the results for a specific source sce- 
nario when ~ 1 /50 of the ORS galaxies more luminous than 
M\ lm = -20.5 is randomly selected as UHECR sources. How- 
ever, the statistical quantities presented in this section are cal- 
culated with not only the statistical error but also the variation 
between different selections of source from our ORS sample. 

The upper panels of figure 5 shows the first and the sec- 
ond harmonics predicted by our source model as a function of 
the cosmic-ray energies for Bsoiar = 1-5/u G, where Z?s iar is me 
strength of the GMF in the vicinity of the Solar system. It is 
noted that we calculate the harmonic amplitudes for the simu- 
lated events within only -10° < S < 80° in order to compare 
with AGASA data. The errorbars represent the statistical fluc- 
tuations due to the finite number of the simulated events, which 
is set to be equal to that observed by the AGASA (775 events 
for E > 10 19 eV). The event selections are performed 20 times. 
The shaded regions represent 1 a total error due to not only 
the statistical error but also the source selections from our ORS 
sample. The random source selections are performed 100 times. 



The region below the histogram is expected values for the statis- 
tical fluctuation of isotropic source distribution with the chance 
probability larger than 10%. 

It is clear that our source model predicts the large-scale 
isotropy fully consistent with that expected by uniform source 
distribution within 1 a total error (statistical one plus source se- 
lection). We have checked that 27 source distributions out of 
100 predict the sufficient large-scale isotropy within 1 a statis- 
tical error. In order to investigate the effects of the GMF on 
the large-scale anisotropy, we also calculate the harmonic am- 
plitude for the case of Bsoiar = 0.0/x G. For Z?s iar = 0.0/1 G, the 
predicted arrival distribution is relatively more isotropic than 
that for Bsoiar = 1-5 At G. We also note that this tendency can 
be seen at lower energies (~ 10 19 eV). Because the deflection 
angle of cosmic rays with such energies by the GMF is about 
~ afew x 10°, the harmonic amplitude of arrival distribution of 
UHECRs can be affected by anisotropy of the event distribu- 
tions which is caused by the events aligned according to the 
order of their energies, reflecting the direction of the GMF. 
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FIG. 6. — Two point correlation function expected for our source model 
for E > 4 X 10 19 (left) and E > 10 19 eV (right). The errorbars represent the 
statistical fluctuations due to the finite number of the simulated events, which 
is set to be equal to that observed by the AGASA within -10° < <5 < 80°. 
The shaded regions represent 1 a total error due to not only the statistical error 
but also the source selections from our ORS sample. The histograms represent 
the AGASA data. However, the AGASA data for E > 10 19 eV are fitted to 
the result of our calculation at larger angle (30° ), since we can not know the 
normalization of the AGASA data with this energy. 

In figure 6, we show two point correlation function predicted 
by our source model for E > 4 x 10 19 (left) and E > 10 19 eV 
(right). It is noted that we calculate two point correlation func- 
tion for the simulated events within only -10° < 5 < 80° in or- 
der to compare with AGASA data. The errorbars represent the 
statistical fluctuations due to the finite number of the simulated 
events, which is set to be equal to that observed by the AGASA 
(775 events for E > 10 19 eV). The shaded regions represent 1 
a total error due to not only the statistical error but also the 
source selections from our ORS sample. The event selections 
and the random source selections are performed 20 times and 
100 times, respectively. The event numbers shown in this fig- 
ures are averaged over all trials of the event selections and the 
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• 10 190 eV<E<10 19/l eV 
10 191 eV<E<10 192 eV 

• 10 192 eV<E<10 194 eV 



10 194 eV<E<10 196 eV 
• 10 19 - 6 eV<E 



arrival directions (earth) 



source directions 



360 c 




FIG. 4. — Arrival directions of protons with E > 10 eV at the earth (left panel) expected for the source model of figure 2 in the galactic coordinate. The events 
are shown by color according to their energies. The event number within —10° < 6 < 80° is set to be equal to the one observed by the AGASA (Takeda et al. 2001, 
775) (The total number of events is ~ 1500). The right panel is the mapping of the sources which actually give rise to events shown in the left panel. Note that this 
mapping differ from the distribution of the sources shown in figure 2. 
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FIG. 5. — Harmonic amplitude predicted by our source model as a function of the cosmic-ray energies. The errorbars represent the statistical fluctuations due to 
the finite number of the simulated events, which is set to be equal to that observed by the AGASA within —10° < S < 80°. The shaded regions represent 1 a total 
error due to not only the statistical error but also the source selections from our ORS sample. The region below the solid line is expected values due to the statistical 
fluctuation of isotropic source distribution with the chance probability larger than 10%. 
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source selections. The histograms represent the AGASA data. 
However, the AGASA data for E > 10 19 eV are fitted to the 
result of our calculation at larger angle (30°), since we can not 
know the normalization of the AGASA data with this energy. 

Clearly visible is that large peak at small angle scale is too 
strong compared with the AGASA observation (Takeda et al. 
2001). We have checked that when extremely nearby sources 
are selected by accident, predicted small-scale anisotropy be- 
comes to be very strong. This is the reason for too large peak at 
small angle scale in figure 6, where the averages and the vari- 
ances are calculated including such source distributions. 

Provided that extremely nearby sources are selected by acci- 
dent, not only the small-scale anisotropy but also the large-scale 
isotropy is inconsistent with the AGASA observation. Accord- 
ingly, we calculate two point correlation function only for the 
source distributions which predict the large-scale isotropy con- 
sistent with uniform source distribution within 1 a statistical 
error. The number of such source distributions is 27 out of all 
the 100 source selections. The result is shown in figure 7. 

Two point correlation function in figure 7 exhibits a structure 
that is similar to that seen in the AGASA data, that is, large 
peak at small angle scale followed by a tail at large angles. For 
E > 10 19 eV, a peak at small angles is somewhat weaker than 
that for E > 4 x 10 19 eV because of the large deflection by the 
GMF. This feature is also seen in AGASA data (Takeda et al. 
2001). From this result, it can be understood that the source 
distributions which predict sufficient large-scale isotropy also 
predict the small-scale anisotropy that is similar to that seen in 
the AGASA data. 
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FIG. 7. — Same as figure 6. But, this is the result only for source distribu- 
tions which predict the large-scale anisotropy consistent with uniform source 
distribution within 1 a statistical error. 

However, we should note that the peak at small angle scale is 
still relatively strong compared with the AGASA. There may be 
two possible explanations for this fact. First, we neglect the ef- 
fects of the extragalactic magnetic field in this study, in order to 
save the CPU time. If we can include this effect by future stud- 
ies, strong correlation at small scale will be reduced because of 
the deflection of UHECRs in the intergalactic space. Second, 
we also neglect the random component of the GMF in order to 



make easy to see the effect of the regular field. This may also 
relax the large peak at small angle scale, provided that there is 
the random component with same level of strength with the reg- 
ular component. These issues are left for future investigations. 

4.3. Future Prospects ofUHECR arrival distribution 

In this subsection, we demonstrate the results of the UHECR 
arrival distribution above E > 10 19 eV with the event num- 
ber expected by future experiments, such as Auger, EUSO and 
OWL. The results for the source model of figure 2 are shown 
in figure 8. The events are shown by color according to their 
energies. It is noted that the expected event rate by the Auger 
experiment is ~ 3000 per year above 10 19 eV. These results are 
extended versions of our previous study (Yoshiguchi, Nagataki, 
& Sato 2003c), where we predict the UHECR arrival distribu- 
tion above 4 x 10 19 eV without modifications by the GMF. 

Remarkable feature is the arrangement of clustered events 
at the directions of nearby sources (see figure 2). The events 
are aligned according to the order of their energies reflecting 
the direction of the GMF. We will be able to obtain some kind 
of information about the GMF and the chemical composition 
of UHECRs. In forthcoming work, we plan studies about new 
statistical quantities which allow us to obtain such invaluable 
information. 

5. SUMMARY AND DISCUSSION 

In this paper, we presented a new method for calculating 
the arrival distribution of UHECRs, which can be applied to 
several source location scenarios, including modifications by 
the GMF. We performed numerical simulations of UHE anti- 
protons, which are injected isotropically at the earth, in the 
Galaxy and recorded the directions of velocities at the earth and 
outside the Galaxy for all of the trajectories. It is noted that we 
regard these anti-protons as PROTONs injected from the out- 
side of the Galaxy toward the earth. We then selected some 
trajectories so that the resultant mapping of the velocity direc- 
tions outside the Galaxy of the selected trajectories corresponds 
to our source location scenario, applying Liouville's theorem. 

There are two points of our improvement over the work of 
Alvarez-Muniz, Engel & Stanev (2002). First, we calculated 
only the trajectories actually reaching the detectors by propa- 
gating anti-protons backwards from Earth, instead of propagat- 
ing protons from the source and selecting those reaching the 
Earth. This helps us to save the CPU time efficiently and makes 
calculations of propagation of cosmic rays even with lower en- 
ergies (<~ 10 19 eV) possible within a reasonable time. Second, 
we considered energy loss processes of UHE protons in the in- 
tergalactic space, which is not taken into account in Alvarez- 
Muniz, Engel & Stanev (2002). This enables us to include the 
effects of difference between resultant energy spectra injected 
at different distances into numerical calculations. We can cal- 
culate the arrival distribution of UHECRs under more realistic 
situation. 

As an application of this method, we calculate the UHECR 
arrival distribution above 10 19 eV for the source model which 
is adopted in our recent study (Paper I) and can explain the cur- 
rent AGASA observation above 4 x 10 19 eV. We found that the 
predicted large-scale anisotropy is fully consistent with uniform 
source distribution, in the same manner as the current AGASA 
data. In order to investigate the effects of the GMF on the 
large-scale anisotropy, we calculated the harmonic amplitude 
for the case of Bsoiar = O.Ofi G. For fisoiar = 0.0 fi G, the pre- 
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FIG. 8. — Arrival directions of protons with E > 10 9 eV at the earth expected for the source model of figure 2 in the galactic coordinate. The events are shown 
by color according to their energies. It is noted that the expected event rate by the Auger experiment is ~ 3000 per year above 10" eV. 



dieted arrival distribution is relatively more isotropic than that 
for flsoiar = 1 .5/it G. This would be due to anisotropy of the event 
distributions which is caused by the events aligned according to 
the order of their energies, reflecting the direction of the GMF. 

It is also found that the calculated two point correlation func- 
tion is similar to that of AGASA data, when we restrict our 
attention to the source distributions which predict sufficient 
isotropic arrival distribution of UHECRs. There may be effects 
of the extragalactic magnetic field and the random component 
of the GMF on the large peak of two point correlation function. 
These issues are left for future studies. 

Finally, we demonstrated the UHECR arrival distribution 
above 10 19 eV with the event number expected by future ex- 
periments in the next few years. The interesting feature of the 
resultant arrival distribution is the events aligned according to 
the order of their energies, reflecting the directions of the galac- 
tic magnetic field. This is also pointed out by Alvarez-Muniz, 
Engel & Stanev (2002). This feature will become clear with in- 
creasing amount of data, and allow us to obtain some kind of in- 
formation about the composition of UHECRs and the GMF. In 
forthcoming work, we plan studies about new statistical quan- 
tities which allow us to obtain such invaluable information. 

In the present work, we calculate the arrival distribution of 



UHECRs for our source location scenario which is adopted in 
our previous study (Paper I). However, it should be mentioned 
that the same results would be obtained if the sources were truly 
drawn at random, as far as the source number density is ~ 10~ 6 
Mpc" 3 . The results such as the events aligned according to the 
order of their energies are independent on our assumption about 
the source distribution. 

In this study, we calculate the harmonic amplitude and two 
point correlation function, which are only published quanti- 
ties on UHECRs observed by the AGASA including lower en- 
ergy ones (~ 10 19 eV). In particular, the AGASA observation 
has published neither existence nor non-existence of the events 
aligned according to the order of their energies. If more de- 
tailed event data with E < 4x 10 19 eV are published, we may 
obtain more strong constraint on our source model, other than 
the source number density, using another statistical quantities. 
This is also one of future study plans. 

The authors acknowledge fruitful discussions with S. Tsub- 
aki. This research was supported in part by Giants-in-Aid 
for Scientific Research provided by the Ministry of Educa- 
tion, Science and Culture of Japan through Research Grant 
No.S14102004 and No. 14079202. 
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